#include "FEMFunction.h"
#include "Equation.h"
#include "BoundaryFunction.h"
#include "FEMSpace.h"
#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi * p[0]) * sin(pi * p[1]);
}
int main()
{
    RectangleDomain* r = new RectangleDomain({{0,0},{1,0},{1,1},{0,1}});
    std::vector<Boundary<2> > B = r->boundary();
    int segmentx = 8;
    int segmenty = segmentx;
    Mesh<2>* m = new Q2Mesh(r,{segmentx,segmenty});
    Element<2>* e = new Quadrilateral_2_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ ->SetRightHandsTermFunction(f);
    BoundaryFunction<2> * bf = new Dirichlet<2>(bc,B);
    BoundaryCondition<2> bfc;
    bfc.add(bf);
    Possion_2D possionproblem(m,e,equ,bfc);
    possionproblem.AssembleStiffMatrix();
    possionproblem.AssembleRightHandsTerm();
    possionproblem.DealWithBoubdaryCondition();
    possionproblem.Solve();
    std::cout << possionproblem.solution() << std::endl;
    FEMFunction<2> v(possionproblem.solution());
    for(int k = 0;k < m->n_element();k++)
    {
        std::vector<int> Idx = m->NodeofEle(k);
        std::vector<Dofs<2>> temp(e->n_Dofs());
        for(int i = 0;i < e->n_Dofs();i++)
            temp[i] = m->DofsofIndex(Idx[i]);
        e->SetDofsList(temp);
        Point<2> p = {0.5,0.5};//the point in reference element
        Dofs<2> p_r = {e->Global_x(0.5,.5),e->Global_y(0.5,.5)};       
        std::cout << bc(*p_r) - v.value(p,e) << std::endl;//计算一个误差
        std::cout << "( " << v.gradient(p,e)[0] << ", " <<  v.gradient(p,e)[1] << " )" <<std::endl;
    }

}